2020地理设计组三等奖：基于深度学习的气象干旱时空预测分析

寒假好！

今天跟大家分享的是2020易智瑞GIS开发竞赛地理设计组三等奖的一个作品：

作品编号：B689（三等奖）

作品名称：基于深度学习的气象干旱时空预测分析

作者单位：华北水利水电大学地球科学与工程学院

小组成员：张棋，丁严，李缙昌

指导老师：许德合

作品视频

作品简介

一、背景与简介从全球范围看，旱灾影响面积最广，造成经济损失最大，被认为是世界上最严重的自然灾害类型之一[1]，而中国干旱地区面积占全国自然灾害影响面积的60%[2]。由于其出现次数多、持续时间长、影响范围大、对农业等经济部门造成的直接损失重，再加上水资源、土地资源及其对社会的潜在影响，成为我国最大的自然灾害之一[3]。因此，量化研究干旱时空变化特征并阐述其形成的机制，对农业生态系统的科学管理和气象灾害预警等方面有重要意义[4]。国际上将干旱分为气象干旱、农业干旱、水文干旱和社会经济干旱四类[5-6]。其中，气象干旱的评估是农业干旱、水文干旱和社会经济干旱监测和预警的基础。通常我们选用便于计算的干旱指标来监测评估干旱发生的强度、持续时间和受灾范围 [7]。目前，气象干旱常用指标主要有标准化降水指数 （Standard Precipitation Index, SPI）、帕默尔干旱指数（Plamer Drought Severity Index, PDSI）、标准化降水蒸散指数（Standard Precipitation Evaporation Index, SPEI）等指数[6]。其中，SPEI适用范围广，易于计算，因此，成为在实际运用当中最广泛并且适用于所有气候的干旱指标[7]。Yao Ning等[8]利用SPEI等干旱指数按照中国7大分区评估了1961-2013年中国大陆干旱的时空演变，并结合SPI等干旱指数揭示了历史的干旱状况。N.S.Abeysingha等人[9]利用SPEI在3个月、6个月和12个月的时间尺度上评估了斯里兰卡的干旱状况，利用曼肯徳尔秩和检验（Mann-Kendall test，M-K）和泰森多边形结合的方式来监测干旱事件发生的趋势及空间分布。李勤[6]等人利用气象干旱和农业干旱结合的综合监测方法对中国1979-2015年不同时间尺度气象干旱和农业干旱进行研究，发现综合监测方法优于单一的气象干旱监测与农业干旱监测。于家瑞[10]等人利用不同时间尺度的SPI值并结合游程理论来分析黑龙省干旱的时空演变特征。综上，对于干旱的监测，目前大多数研究均从时间和空间两个纬度分别进行探究，而未从总体的时空尺度上去分析干旱的特征与变化。利用时空数据模型探索气象干旱的分布格局、形成过程和影响机制，有重要的实践与现实意义[11]。从时空数据方面；本项目利用时空立方体模型对中国6个指数（PDSI、PMDI、PHDI、Z-index,、sc-PDSI、SPEI）从时空尺度上进行可视化展示与分析。对中国613个站点创建泰森多边形，并结合时空尺度的K-means聚类方法对近40年来613个区域进行聚类，利用伪F统计量评估三个时间尺度（SPI3, SPI6和SPI12）的最佳聚类数。对于历史旱情严重年份则利用经验贝叶斯克里金插值（Empirical Bayesian Kriging, EBK）法对其进行空间分布预测研究。从时间序列方面，结合自回归移动平均模型（Autoregressive Integrated Moving Average, ARIMA）和长短时记忆神经网络（Long Short-Term Memory，LSTM），以及二者的组合模型（ARIMA-LSTM）对多尺度指数进行预测。

二、研究区概况与数据

中国地形由西至东呈下降趋势，且地貌类型多样化，由平原、山谷、丘陵、水系、水体、高原、沙漠和冰川盆地等组成[8]。气候类型可分为季风气候、温带大陆性气候、高寒气候，温度带由南至北可划为热带、亚热带、暖温带、中温带、寒带[12]。本文选取1980-2019年中国613个连续监测的气象站逐月降水量数据进行计算，来源于国家气象信息中心提供的中国地面气候资料月值数据集（http://data.cma.cn/）。考虑到中国多年的平均气温条件，中国被划分7个自然分区（见图1），分别为西北荒漠地区（图1中标注为“1”），内蒙古草原地区（图1中标注为“2”），东北湿润半湿润温带地区（图1中标注为“3”），华北湿润半湿润暖温带地区（图1中标注为“4”），华中华南湿润亚热带地区（图1中标注为“5”），青藏高原（图1中标注为“6”），华南湿润热带地区（图1中标注为“7”）。

图1  中国7大地理分区 与612个气象站点分布

三、设计思想

3.1  6个干旱指数（PDSI、PMDI、PHDI、Z-index、sc-PDSI、SPEI）

由于中国各地旱情均有季节区分并且在旱情严重地区均有连旱情况发生，而大部分农作物生长季在3月到8月，因此本文分别计算了1980-2019年的1, 3, 6, 9, 12和24个月的SPEI值以及PDSI、PMDI、PHDI、Z-index、sc-PDSI的值，并通过国家标准气象干旱等级（GB/T20481-2006）规定的干旱分级标准（表1）来表征干旱情况[3]。

表1 标准化降水指数干旱分级

Tab.1 Drought classification based on SPI

等级	类型	指数范围
1	无旱	指数≥-0.5
2	轻旱	-1.0≤指数<-0.5
3	中旱	-1.5≤指数<-1.0
4	重旱	-2.0≤指数<-1.5
5	特旱	指数≤-2.0

3.2  时空模式挖掘

结合时空分析思维模式来探究气象干旱在中国的时空分布特征、时空演化过程、时空聚类分析和时空热点分析，能够为有关部门防旱抗旱提供科学依据。以下研究方法均从时空尺度出发并结合SPI对中国近40年气象干旱时空特征进行研究与分析。

3.2.1  时空立方体     

时空立方体模型是通过将样本点聚合到时空条柱的方法，形成NetCDF（Network Common Data Form）,即网络通用数据格式的时空数据结构中[18]。通过创建时空立方体，能以时间序列分析、集成空间和时间模式分析等形式，对时空数据进行可视化与分析[18]（见图2）。图2中x和y轴代表该时间段的空间位置，z轴代表时间，且最下面的为起始时间，最上面的为最近时间，每个立方体均由该时间对应的属性值组成，数值的大小可以通过设置不同颜色来区分。本文通过将气象站点数据计算得到的SPI值进行聚合，每年的SPI值均聚合到一个立方体中，而多个立方体按照时序排列（垂直排列）并聚合后得到时空柱，即每个站点形成一个时空柱。因此通过聚合不同时间尺度的SPI值，可得到中国612个气象站点监测到的不同类型的旱情。

图2  时空立方体模型 

3.2.2  时空序列聚类    

在时空立方体的基础上，基于时间序列的相似特征，将时空立方体中的时间序列集合进行划分，也就是K-means算法的时空体现，其中每个聚类的成员具有的时间序列特征均相似。可以基于三个条件聚集时间序列：具有相似的时间值，趋于同时增加和减少以及具有相似特征的重复模式[19]。本文对1980-2019年SPI3、SPI6和SPI12结合时空立方体模型进行时空聚类，首先对中国612个气象站点创建泰森多边形，同时结合K近邻算法和K-means算法对612个泰森多边形进行聚类，聚类数目利用伪F统计量来选取，伪F统计量是一个反应聚类间方差和距离内方差的比率，即反应组内相似性和组间差异的比率。其值越大，表示该聚类数越有效。公式如下：

其中：

n=要素数，=聚类i中的要素数，=类（聚类）数，=用于聚类要素的变量数，=i聚类中j要素的k变量值，=k变量的平均值，=聚类i中k变量的平均值。

3.2.3  时空热点分析    

时空热点分析可识别数据的时空趋势，探测某一特征在时空尺度的热点或冷点，在本文中，热点代表该地区往年未有旱情或旱情不明显，而近年发生旱情或相比往年旱情严重的情况，而冷点则相反。通过特定的邻域距离和临域时间步长参数来计算每个立方体条柱的Getis-Ord Gi*统计量[20-21]。结合M-K检验法对时空尺度的热点分析结果进行趋势检验，共有17种，分别为新增的、连续的、加强的、持续的、逐渐减少的、分散的、振荡的以及历史的热点和冷点[21]，还有一个是无显著特征。由于研究方向不同，因此本文共检验出6种，分别为：新兴热点、分散热点、振荡热点、新兴冷点、振荡冷点和无显著特征。

3.2.4  经验贝叶斯克里金插值法（EBK）     

EBK是一种地统计插值方法，专为克服经典克里金法更加困难的理论和实践局限而开发的，它可以用构造子集和模拟的方法自动计算这些参数，而不需要像其它克里金插值法那样手动调整参数，并且可以通过估计基础半变异函数来说明引入的误差，因此EBK法与其他克里金相比预测精度更高[22]。

3.3  时间序列预测

3.3.1  ARIMA模型

ARIMA分为自回归模型(Autoregressive model，AR)、滑动平均模型(Moving average model, MA)以及自回归移动平均模型(Autoregressive Moving Average Model, ARMA)，是传统的时间序列预测模型。其建模流程是首先判断模型平稳程度，其次利用差分法对非平稳时间序列进行平稳化处理，然后选取AR(p),MA(q)对模型进行定阶，差分次数记为d，ARIMA(p,d,q)模型就是经过了d阶差分后的ARMA（p，q）模型。如下式所示： 

（4）中θ(L)和ϕ(L)分别为：

ARIMA(p，d，q)模型的一般式：

上式中，d为差分次数，为自回归系数，p为自回归阶数，为移动平均系数；为白噪声序列（服从0均值、正态分布且相互独立的白噪声序列）。

其中，p，q阶数采用赤池信息准则（Akaike Information Criterion, AIC）和贝叶斯信息准则（Bayesian Information Criterion, BIC）来确定，当样本数N固定时，选择AIC和BIC最小值来确定p，q。公式如下：

3.3.2  LSTM模型

LSTM是RNN的一种拓展，结构见图2。LSTM与RNN最大的不同在于隐藏模块中添加了三个门层，分别是忘记门（Forget gate）、输入门（Input gate）和输出门（Output gate）。忘记门决定什么样的信息保留，什么样的信息遗忘；输入门用来更新细胞状态，决定什么值我们需要更新；输出门是一个过滤器，基于我们的细胞状态来选择哪个部分用来输出。其中LSTM隐藏层结构如图3所示，其计算方法可以表示为[27]：

式中：i、f、C、o分别表示为输入门、遗忘门、细胞状态、输出门；W和b分别为对应的权重系数矩阵和偏置项；σ和tanh分别为sigmoid和双曲正切激活函数。

图 3 LSTM模型结构图

图 4 LSTM隐藏层细胞结构

本文利用LSTM对多时间尺度SPI序列进行拟合，由于LSTM的隐含层数量决定了模型拟合能力，为了防止过拟合，本文提出了“早停法”，即当损失函数不再下降时停止训练。损失函数定义如下：

式中：在时间i处的观测值，是在时间i处的预测值。

建模流程如下：

·输入数据的预处理

通常，神经网络模型处理的数据是归一化数据，范围在[-1,1]。因为归一化可以使学习率不必再根据数据范围做调整，从而提高模型的训练速度。

·参数率定与训练

本文基于Python3.7平台来搭建LSTM网络，采用基于时间的反向传播算法进行训练。batch_size设置为1，其中1代表在每个样本之后更新权重，该过程称为随机梯度下降（Stochastic Gradient Descent, SGD）。激活函数通常有sigmoid、tanh和relu，由于sigmoid在反向传播时很容易出现梯度消失情况，而tanh的SGD收敛速度太慢，因此本文选用relu做激活函数。由于本文选用“早停法”来防止训练过拟合，该方法选取均方误差（Mean Squared Error, MSE）作为损失函数，损失函数越小，说明预测模型精度越高，随着迭代次数不断增加时MSE值逐渐下降，直到某一时刻MSE值上升时说明模型过拟合，因此早停法在模型MSE值上升前一刻停止训练以保证预测模型的精度达到最高。迭代次数本文设置为300次，以保证MSE能降到最低。对于LSTM网络而言，神经元隐藏层数量决定着网络训练速度与预测精度，隐藏层数量太少，导致网络无法训练或效果很差，隐藏层数量过大又导致网络训练速度慢或者出现过拟合现象。本文采用黄金分割法来选择隐藏层神经元数量。黄金分割法是指在区间[a，b]中找到理想的隐含层节点个数，然后根据黄金分割原理扩展搜索区间，即得到区间[b，c]（b=0.619*（c-a）+a）。

·网络输出与评价验证

由于数据在预处理时采用归一化处理，因此对LSTM模型残差预测完毕后，得到的并不是最终结果，还需将该时间序列数据进行反归一化处理，即得到模型实际输出。经过交叉验证得到80%（1951.1—2004.8）数据作为训练集，20%（2004.9—2017.12）数据作为测试集。

3.3.3  ARIMA-LSTM模型

由于ARIMA模型和LSTM模型在线性和非线性预测中各有优点，因此本文分别采用ARIMA与LSTM模型的优点建立组合模型ARIMA-LSTM,假设时间序列可视为线性自相关部分与非线性残差两部分的组合，利用ARIMA模型对SPI值进行预测，将结果与实际值相减得到残差，将残差记为非线性部分，将残差带入LSTM模型进行预测，最后把两者预测结果相加得到组合结果，即：

四、主要功能本研究的主要功能是结合线性和非线性模型对多尺度干旱指数进行组合预测，并利用改进的克里金插值法EBK 和时空模式挖掘方法对预测结果进行可视化分析与展示。能够为我国气象干旱、农业生产，防灾减灾提供理论基础与科学依据。

五、特点与创新本项目的特点与创新之处在于：（1）针对现有时间序列预测方法要么是线性模型，要么是非线性模型，而线性和非线性模型在不同时间尺度上的预测效果各有利弊，提出将线性模型和非线性模型结合的方法：本项目通过对不同时间尺度的SPI 序列进行线性模型和非线性模型的组合预测，解决了大部分时间序列预测问题时只能“二选一”的难题，同时本文在非线性模型上选用了具有自学习能力的LSTM 来替代传统的ANN 等非线性模型，解决了神经网络在训练过程中极大可能会产生过拟合（overfit）和欠拟合（underfit）的现象，从而提高非线性模型的预测精度。（2）针对现有干旱的预测分析中，绝大多数研究仅对时间序列进行预测，并未考虑其时空关系，因此本项目采用深度学习中基于神经网络的空间组合来捕捉各站点空间相关性，从而对时空序列进行预测。（3）针对现有干旱的时空分析中，绝大多数研究仅对个别年份进行二维空间分析与展示，不能从三维角度分析多年时空分布的情况，在此基础上提出EBK 与时空立方体结合的方法：本项目通过EBK 与时空立方体模型相结合对预测结果进行时空分布展示与分析。

相关
文章

精彩推荐
分类	内容	适宜对象
课程篇	ArcGIS Pro制图入门系列课程	在读生 易智瑞开发竞赛选手
课程篇	Web GIS课程全链接	在读生 GIS应用开发组选手
课程篇	Portal for ArcGIS入门	有GIS开发基础 GIS应用开发组选手
考研系列	写在考研成绩之时（1）选导师	考研、复试同学
考研系列	写在考研成绩之时（2）联系导师	考研、复试同学
考研系列	写在考研成绩之时（3）复试	考研、复试同学
考研系列	来，一起扒一下地信专业考研数据	准备考研的同学
学习篇	这是一份遥感学习笔记：遥感数据源&免费数据获取途径	在读生 遥感应用组同学
学习篇	人生苦短，你需要Python（0基础患者入门篇）	Python入门 地理设计组同学
学习篇	人生苦短系列2——来自ArcGIS Python的拯救	Python入门 地理设计组同学
学习篇	带你做一遍“垃圾场选址”这个case（一）	地理设计组同学
学习篇	带你做一遍“垃圾场选址”这个case（二）	地理设计组同学
学习篇	带你做一遍“垃圾场选址”这个case（三）	地理设计组同学
竞赛篇	用地图，讲故事，2020开发竞赛，我们认真的！	地图故事组同学
竞赛篇	ArcGIS Online中国 上的『俄罗斯三城游记』	地图故事组同学
竞赛篇	如何做出好作品？开发竞赛导师宝藏箴言（完整版）	参赛选手
竞赛篇	北京大学刘瑜教授：能解决实际问题才是真的GIS人才	参赛选手
竞赛篇	“C”位出道，开发竞赛助你成为下一个开挂的GISer	参赛选手

关注我们

公众号ID：M姐实验室

关注GIS教育的一切